FA=function(Y,K=2,sigma,epsilon = 1e-3,max_iter=200){
set.seed(601)
n=dim(Y)[1]
p=dim(Y)[2]
lambda=matrix(c(rep(1,4),numeric(3),numeric(3),rep(1,4)),nrow=7,ncol=2)+matrix(rnorm(p*K,0,sigma),nrow=p,ncol=K)
phi=diag(rnorm(p,0,1))
loss=c()
original=matrix(c(rep(1,4),numeric(3),numeric(3),rep(1,4)),nrow=7,ncol=2)%*%t(matrix(c(rep(1,4),numeric(3),numeric(3),rep(1,4)),nrow=7,ncol=2))
new=lambda%*%t(lambda)
diff=norm(new-original,type="F")
t=0
while(diff>epsilon){
t=t+1
Varx=diag(1,K)-t(lambda)%*%solve(lambda%*%t(lambda)+phi)%*%lambda
Ex=Y%*%solve(lambda%*%t(lambda)+phi)%*%lambda
Extx=n*Varx+t(Ex)%*%Ex
lambda_new=t(Y)%*%Ex%*%solve(Extx)
lambda=lambda_new
phi_new=diag(diag(1/n*(t(Y)%*%Y-2*t(Y)%*%Ex%*%t(lambda)+lambda%*%Extx%*%t(lambda))))
phi=phi_new
Sigma_y=lambda%*%t(lambda)+phi+epsilon * diag(p)
loss<-c(loss,(-n*p/2)*log(2*pi)-(n/2)*log(abs(det(Sigma_y)))-0.5*sum(diag(solve(Sigma_y)%*%t(Y)%*%Y)))
new=lambda%*%t(lambda)
diff=norm(new-original,type="F")
if(t>max_iter){
break
}
}
#print(paste("Times of iteration=",i))
return(list(lambda=lambda,phi=phi,loss=loss))
}
n=100
generator_pca=function(n=100){
Lambda=matrix(c(rep(1,4),numeric(3),numeric(3),rep(1,4)),nrow=7,ncol=2)
X=mvrnorm(n,numeric(2),diag(1,2))
W=mvrnorm(n,numeric(7),diag(0.4,7))
Y=t(Lambda%*%t(X)+t(W))
return(Y)
}
set.seed(601)
Y=generator_pca(100)
A1=FA(Y,2,1,1e-3,70)
#set.seed(601)
A2=FA(Y,2,3,1e-3,70)
#set.seed(601)
A3=FA(Y,2,5,1e-3,70)
#set.seed(601)
A4=FA(Y,2,7,1e-3,70)
A5=FA(Y,2,9,1e-3,70)
A6=FA(Y,2,11,1e-3,70)
df1=NULL
df1=data.frame(iteration=c(1:length(A1$loss)),loss=A1$loss)
df2=data.frame(iteration=c(1:length(A2$loss)),loss=A2$loss)
df3=data.frame(iteration=c(1:length(A3$loss)),loss=A3$loss)
df4=data.frame(iteration=c(1:length(A4$loss)),loss=A4$loss)
df5=data.frame(iteration=c(1:length(A5$loss)),loss=A5$loss)
df6=data.frame(iteration=c(1:length(A6$loss)),loss=A6$loss)
plot_ly(data=df1,x=~iteration,y=~loss,type="scatter",mode="lines+markers",name = "Sigma=1")%>%
add_trace(x = df1[,1], y = df2[,2],
type = "scatter", mode = "lines+markers",name = "Sigma=3")%>%
add_trace(x = df1[,1], y = df3[,2],
type = "scatter", mode = "lines+markers",name = "Sigma=5")%>%
add_trace(x = df1[,1], y = df4[,2],
type = "scatter", mode = "lines+markers",name = "Sigma=7")%>%
add_trace(x = df1[,1], y = df5[,2],
type = "scatter", mode = "lines+markers",name = "Sigma=9")%>%
add_trace(x = df1[,1], y = df6[,2],
type = "scatter", mode = "lines+markers",name = "Sigma=11")